clear ;
global gamma bet Bx dx dt 

dx = 0.008 ;
dt = 0.000025 ; 
gamma = 5/3;
bet = 2;
tstep = 200 ;
X = (-1:dx:1) ;
Bx = 5;

rho = zeros(length(X),tstep) ;
p = zeros(length(X),tstep) ;
vx = zeros(length(X),tstep) ; 
vy = zeros(length(X),tstep) ; 
vz = zeros(length(X),tstep) ; 
By = zeros(length(X),tstep) ; 
Bz = zeros(length(X),tstep) ; 
Q = zeros (7,3); 
M = zeros(7,2) ; 

%initlization ==========

for i = 1:length(X) 
    if i <= (length(X)-1)/2+1
        rho(i,1) = 3.108 ;
        p(i,1) = 1.4336 ;
        vx(i,1) = 0 ;
        vy(i,1) = 0.2633 ;
        vz(i,1) = 0.2633 ;
        By(i,1) = 0.1 ;
        Bz(i,1) = 0.1 ;
    else
        rho(i,1) = 1.0 ;
        p(i,1) = 0.1 ;
        vx(i,1) = -0.9225 ;
        vy(i,1) = 0.0 ;
        vz(i,1) = 0.0 ;
        By(i,1) = 1.0 ;
        Bz(i,1) = 1.0 ; 
    end 
end
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 
for n = 1 : tstep 
    n
    rho(2,n)=rho(3,n) ;
    p(2,n)=p(3,n) ;
    vx(2,n)=vx(3,n) ; 
    vy(2,n)=vy(3,n) ;
    vz(2,n)=vz(3,n) ;
    By(2,n)=By(3,n) ;
    Bz(2,n)=Bz(3,n) ;
    
    rho(1,n)=rho(2,n) ;
    rho(length(X),n) = rho(length(X)-1,n) ; 
    p(1,n)=p(2,n) ;
    p(length(X),n) = p(length(X)-1,n) ; 
    vx(1,n)=vx(2,n) ;
    vx(length(X),n) = vx(length(X)-1,n) ; 
    vy(1,n)=vy(2,n) ;
    vy(length(X),n) = vy(length(X)-1,n) ; 
    vz(1,n)=vz(2,n) ;
    vz(length(X),n) = vz(length(X)-1,n) ; 
    By(1,n)=By(2,n) ;
    By(length(X),n) = By(length(X)-1,n) ; 
    Bz(1,n)=Bz(2,n) ;
    Bz(length(X),n) = Bz(length(X)-1,n) ; 
    
    for m = 3 : length(X)-1
        Q(1,1) = rho(m-1,n) ; 
        Q(1,2) = rho(m,n) ; 
        Q(1,3) = rho(m+1,n) ;
        
        Q(2,1) = p(m-1,n) ; 
        Q(2,2) = p(m,n) ; 
        Q(2,3) = p(m+1,n) ;
        
        Q(3,1) = vx(m-1,n) ; 
        Q(3,2) = vx(m,n) ; 
        Q(3,3) = vx(m+1,n) ;
        
        Q(4,1) = vy(m-1,n) ; 
        Q(4,2) = vy(m,n) ; 
        Q(4,3) = vy(m+1,n) ;
        
        Q(5,1) = vz(m-1,n) ; 
        Q(5,2) = vz(m,n) ; 
        Q(5,3) = vz(m+1,n) ;
        
        Q(6,1) = By(m-1,n) ; 
        Q(6,2) = By(m,n) ; 
        Q(6,3) = By(m+1,n) ;
        
        Q(7,1) = Bz(m-1,n) ; 
        Q(7,2) = Bz(m,n) ; 
        Q(7,3) = Bz(m+1,n) ;
        
        M(:,1) =( Q(:,1) + Q(:,2) ) / 2;
        M(:,2) =( Q(:,3) + Q(:,2) ) / 2;
        
        flux = rightmat(Q(1,2),Q(2,2),Q(3,2),Q(4,2),Q(5,2),Q(6,2),Q(7,2)) * diag(abs(eigenvalue(Q(1,2),Q(2,2),Q(3,2),Q(4,2),Q(5,2),Q(6,2),Q(7,2)))) * ...
            leftmat(Q(1,2),Q(2,2),Q(3,2),Q(4,2),Q(5,2),Q(6,2),Q(7,2)) * Q(:,2) + negflux(Q(1,3),Q(2,3),Q(3,3),Q(4,3),Q(5,3),Q(6,3),Q(7,3)) - posflux(Q(1,1),Q(2,1),Q(3,1),Q(4,1),Q(5,1),Q(6,1),Q(7,1)) - ...
            secondflux(M(1,1),M(2,1),M(3,1),M(4,1),M(5,1),M(6,1),M(7,1)) * lim(rho,p,vx,vy,vz,By,Bz,m)+ ...
            secondflux(M(1,2),M(2,2),M(3,2),M(4,2),M(5,2),M(6,2),M(7,2)) * lim(rho,p,vx,vy,vz,By,Bz,m+1) ;
        
        
        rho(m,n+1) = rho(m,n) -dt/dx* flux(1) ;
        p(m,n+1) = p(m,n) -dt/dx* flux(2) ;
        vx(m,n+1) = vx(m,n) -dt/dx* flux(3) ;
        vy(m,n+1) = vy(m,n) -dt/dx* flux(4) ;
        vz(m,n+1) = vz(m,n) -dt/dx* flux(5) ;
        By(m,n+1) = By(m,n) -dt/dx* flux(6) ;
        Bz(m,n+1) = Bz(m,n) -dt/dx* flux(7) ;
        
    end 
   
end 


plot(X,rho(:,tstep))



function param = getparam(rho,p,vx,vy,vz,By,Bz)
global bet gamma Bx
%order as:
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 
    param(1) = sqrt (bet * gamma * p /rho) ;
    param(2) = abs(Bx)/sqrt(rho) ;
    param(3) = sqrt (By^2 + Bz^2) / sqrt (rho) ;
    param(4) = sqrt (param(3)^2 + param(2)^2) ; 
    a = param(1) ;
    bx = param(2) ;
    b = param(4) ; 
    param(5) = (0.5*(a^2 + b^2 - sqrt ((a^2 + b^2 )^2 - 4*a^2*bx^2)))^0.5 ; 
    param(6) = (0.5*(a^2 + b^2 + sqrt ((a^2 + b^2 )^2 - 4*a^2*bx^2)))^0.5 ; 
    param(7) = By/sqrt(By^2+Bz^2) ; 
    param(8) = Bz/sqrt(By^2+Bz^2) ; 
    cs = param(5) ;
    cf = param(6) ; 

    param(9) = sqrt ((cf^2-a^2)/(cf^2-cs^2)) ;
    param(10) = sqrt ((a^2-cs^2)/(cf^2-cs^2)) ;
    param(11) = sqrt(vx^2 + vy^2 + vz^2); 
   
end 

function lambda = eigenvalue(rho,p,vx,vy,vz,By,Bz)
%order of eigenvalues vx vxcs vxbx vxcf corresponding to rho p vx vy vz by
%bz 
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 
global bet gamma Bx
    para = getparam(rho,p,vx,vy,vz,By,Bz) ;
    cs = para(5) ;
    cf = para(6) ; 
    bx = para(2) ;
   
    lambda(1) = vx ;
    lambda(2) = vx - cs ; 
    lambda(3) = vx + cs ; 
    lambda(4) = vx + bx ;
    lambda(5) = vx - bx ; 
    lambda(6) = vx - cf ; 
    lambda(7) = vx + cf ;
end 

function R = rightmat(rho,p,vx,vy,vz,By,Bz)
global Bx bet gamma
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 

    para = getparam(rho,p,vx,vy,vz,By,Bz);
    a = para(1) ;
    bx = para(2) ;
    bp = para(3) ; 
    b = para(4) ;
    cs = para(5) ;
    cf = para(6) ; 
    epsy = para(7) ;
    epsz = para(8) ;
    as = para(9) ; 
    af = para(10) ; 
    v = para(11) ;
   
    
    R(:,1) = [1;0;0;0;0;0;0] ;
    
    R(:,2) = [0;bet*as/a^2;-rho*as*cs/a^2;-rho*af*cf*epsy*sign(Bx)/a^2;-rho*af*cf*epsz*sign(Bx)/a^2;...
        -sqrt(rho)*af*epsy/a;-sqrt(rho)*af*epsz/a];
    
    R(:,3) = [0;bet*as/a^2;rho*as*cs/a^2;rho*af*cf*epsy*sign(Bx)/a^2;rho*af*cf*epsz*sign(Bx)/a^2;...
        -sqrt(rho)*af*epsy/a;-sqrt(rho)*af*epsz/a];
    
    R(:,4) = [0;0;0;epsz/sqrt(2);-epsy/sqrt(2);-epsz*sign(Bx)*sqrt(rho/2);epsy*sign(Bx)*sqrt(rho/2)];
    
    R(:,5) = [0;0;0;epsz/sqrt(2);-epsy/sqrt(2);epsz*sign(Bx)*sqrt(rho/2);-epsy*sign(Bx)*sqrt(rho/2)];
    
    R(:,6) = [af/2;a^2*af/2/bet;-af*cf/2/rho;as*cs*epsy*sign(Bx)/2/rho;as*cs*epsz*sign(Bx)/2/rho;...
        a*as*epsy/2/sqrt(rho); a*as*epsz/2/sqrt(rho)]; 
    
    R(:,7) = [af/2;a^2*af/2/bet;af*cf/2/rho;-as*cs*epsy*sign(Bx)/2/rho;-as*cs*epsz*sign(Bx)/2/rho;...
        a*as*epsy/2/sqrt(rho); a*as*epsz/2/sqrt(rho)];     
end 
   
function L = leftmat(rho,p,vx,vy,vz,By,Bz)
global Bx bet gamma
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 

    para = getparam(rho,p,vx,vy,vz,By,Bz) ;
    a = para(1) ;
    bx = para(2) ;
    bp = para(3) ; 
    b = para(4) ;
    cs = para(5) ;
    cf = para(6) ; 
    epsy = para(7) ;
    epsz = para(8) ;
    as = para(9) ; 
    af = para(10) ; 
    v = para(11) ;
   
    L(1,:) = [1,-bet/a^2,0,0,0,0,0]; 
    
    L(2,:) = [0,bet*as/a^2 , -rho*as*cs/a^2 , -rho*af*cf*epsy*sign(Bx)/a^2 , -rho*af*cf*epsz*sign(Bx)/a^2,...
        -sqrt(rho)*af*epsy/a, -sqrt(rho)*af*epsz/a];
    
    L(3,:) = [0,bet*as/a^2 , rho*as*cs/a^2 , rho*af*cf*epsy*sign(Bx)/a^2 , rho*af*cf*epsz*sign(Bx)/a^2,...
        -sqrt(rho)*af*epsy/a, -sqrt(rho)*af*epsz/a];
    
    L(4,:) = [0,0,0,epsz/sqrt(2), -epsy/sqrt(2) , epsz*sign(Bx)*sqrt(rho*2) , -epsy*sign(Bx)*sqrt(rho*2)];
    
    L(5,:) = [0,0,0,epsz/sqrt(2), -epsy/sqrt(2) , -epsz*sign(Bx)*sqrt(rho*2) , epsy*sign(Bx)*sqrt(rho*2)];
    
    L(6,:) =  [0,bet*af/a^2,-rho*af*cf/a^2,rho*as*cs*epsy*sign(Bx)/a^2,rho*as*cs*epsz*sign(Bx)/a^2,sqrt(rho)*as*epsy/a,sqrt(rho)*as*epsz/a];
    
    L(7,:) =  [0,bet*af/a^2,rho*af*cf/a^2,-rho*as*cs*epsy*sign(Bx)/a^2,-rho*as*cs*epsz*sign(Bx)/a^2,sqrt(rho)*as*epsy/a,sqrt(rho)*as*epsz/a];
end 

function pflux = posflux (rho,p,vx,vy,vz,By,Bz)
global Bx bet gamma
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 
    lambda = eigenvalue(rho,p,vx,vy,vz,By,Bz); 
    R = rightmat(rho,p,vx,vy,vz,By,Bz) ; 
    L = leftmat (rho,p,vx,vy,vz,By,Bz);
    for i = 1:7
        if lambda(i) < 0
            lambda(i) = 0;
        end 
    end 
    poseig = diag(lambda) ;
    uvec = [rho;p;vx;vy;vz;By;Bz];
    pflux = R * poseig * L * uvec ;
end 

function nflux = negflux(rho,p,vx,vy,vz,By,Bz)
global Bx bet gamma
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 
    lambda = eigenvalue(rho,p,vx,vy,vz,By,Bz) ; 
    R = rightmat(rho,p,vx,vy,vz,By,Bz) ; 
    L = leftmat (rho,p,vx,vy,vz,By,Bz) ;
    
    for i = 1:7
        if lambda(i) > 0
            lambda(i) = 0;
        end 
    end 
    negeig = diag(lambda) ;
    uvec = [rho;p;vx;vy;vz;By;Bz] ;
    nflux = R * negeig * L * uvec ;
end 

function secflux = secondflux(rho,p,vx,vy,vz,By,Bz)
global Bx bet gamma dx dt 
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 
    lambda = eigenvalue(rho,p,vx,vy,vz,By,Bz); 
    R = rightmat(rho,p,vx,vy,vz,By,Bz) ; 
    L = leftmat (rho,p,vx,vy,vz,By,Bz);
    
    absA = R * diag(abs(lambda)) * L;
    secflux = 0.5 * absA * (eye(7) - dt / dx * absA); 
end 

function limiter = lim(rho,p,vx,vy,vz,By,Bz,m)
global Bx bet gamma dx dt 
% 1a 2bx 3bp 4b 5cs 6cf 7epsy 8epsz 9as 10af 11v 
    L = leftmat(rho(m),p(m),vx(m),vy(m),vz(m),By(m),Bz(m)) ; 
    L = (L + leftmat(rho(m-1),p(m-1),vx(m-1),vy(m-1),vz(m-1),By(m-1),Bz(m-1))) / 2 ;
    R = rightmat(rho(m),p(m),vx(m),vy(m),vz(m),By(m),Bz(m)) ; 
    R = (R + rightmat(rho(m-1),p(m-1),vx(m-1),vy(m-1),vz(m-1),By(m-1),Bz(m-1))) / 2 ;
    u = [rho(m) - rho(m-1) ; p(m)-p(m-1) ; vx(m)-vx(m-1) ; vy(m) - vy(m-1) ; ... 
        vz(m)-vz(m-1) ; By(m) - By(m-1) ; Bz(m) - Bz(m-1)];
    lambda =  eigenvalue(rho(m),p(m),vx(m),vy(m),vz(m),By(m),Bz(m)) ; 
    lambda = (lambda + eigenvalue(rho(m-1),p(m-1),vx(m-1),vy(m-1),vz(m-1),By(m-1),Bz(m-1))) / 2 ;
    
    alpharray1 = L * u ;
    
    m = m - 1;
    L = leftmat(rho(m),p(m),vx(m),vy(m),vz(m),By(m),Bz(m)) ; 
    L = (L + leftmat(rho(m-1),p(m-1),vx(m-1),vy(m-1),vz(m-1),By(m-1),Bz(m-1))) / 2 ;
    u = [rho(m) - rho(m-1) ; p(m)-p(m-1) ; vx(m)-vx(m-1) ; vy(m) - vy(m-1) ; ... 
        vz(m)-vz(m-1) ; By(m) - By(m-1) ; Bz(m) - Bz(m-1)];

    alpharray0 = L * u ;
    
    m = m + 2;
    L = leftmat(rho(m),p(m),vx(m),vy(m),vz(m),By(m),Bz(m)) ; 
    L = (L + leftmat(rho(m-1),p(m-1),vx(m-1),vy(m-1),vz(m-1),By(m-1),Bz(m-1))) / 2 ;
    u = [rho(m) - rho(m-1) ; p(m)-p(m-1) ; vx(m)-vx(m-1) ; vy(m) - vy(m-1) ; ... 
        vz(m)-vz(m-1) ; By(m) - By(m-1) ; Bz(m) - Bz(m-1)];

    alpharray2 = L * u ;
    limiter = zeros (7,1) ;

    for i = 1 : 7 
        if alpharray1(i) == 0 
            phi = 0 ;
        else 
            if lambda(i) > 0 
                theta = alpharray0(i) / alpharray1(i) ; 
                
            else
                theta = alpharray2(i) / alpharray1(i) ;
            end 

            phi = (theta + abs(theta)) / (1+abs(theta)) ;
        end 
        limiter(:) = limiter(:) + phi * alpharray1(i) * R (:,i) ;
        
    end 

end 
